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I. ESTRODUCTION 



Numerical techniques for solving electromagnetic scattering problems using integral 
equations and the method of moments (MM) are well known. The physical problem, 
specified by Maxwell’s equations and boundary conditions, is reduced to an integro- 
differential equation over finite domains, and solved using a procedure referred to as the 
method of moments [Ref 1]. The unknowns (usually currents) are represented 
by a series of basis functions with unknown expansion coefficients. The MM process 
generates a set of linear equations that must be solved simultaneously. Until recently, 
these techniques have been limited to small (1 to 10 wavelength) geometries because of 
computer run time and memory constraints. With the development of faster computers 
with more memory, the MM technique has increasing application to larger geometries. 
However, computer memory and run time can still be inadequate to solve many 
important antenna and scattering problems. Numerically efficient solutions require less 
computer memory and/or less computer run time. Therefore, any increase in the 
efficiency of a MM solution is of great practical interest. 

The usual MM approach to modeling a thin wire of arbitrary shape is to 
specify a series of points, with piecewise linear segments between the points to 
approximate the wire. The current is represented by one or more basis functions, each 
with constant phase, over a piecewise linear segment. Typically, the size of the segments 
is set by how accurately the current or scattered field needs to be determined. For 
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convergence of the current, segment lengths of 0.05 X to 0.1 X are generally required. 

A second factor that influences the segment size is the radius of curvature of the 
wire. Tightly curved wires require smaller segments to reproduce the wire shape 
accurately. When the radius of curvature is larger than a wavelength, the first case sets 
the segment size; when the radius of curvature is much less than a wavelength, the 
second case dictates the segment size (Figure 1). 

All generally available MM codes based on the method of subdomains use 
the first approach. A natural question arises: Does dividing the wire into curved segments 
that conform to its shape, with arclengths restricted only by the maximum current 
variation rule, yield a solution that converges with fewer subsections? If the answer is 
yes, is the improvement in convergence worth the greater complexity and coding effort? 
To resolve this issue, the following approach is taken: 





Figure 1. Left: Piecewise Linear Segments, Right: Curved Segments 
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1 . 



Formulate the solution for linearly and circularly polarized plane wave incidence 
for a geometrically simple shape such as a loop. 

2. Computer code the solution in FORTRAN. 

3. Validate the solution using other MM solutions and measured data. 

4. Study the convergence of the solution with respect to loop parameters and 
compare its performance to a method using linear subsections. 

5. Study the effect of changes in program structure on its computational efficiency. 
Chapter II discusses the derivation and the MM solution of the electric field integral 

equation for a thin wire loop. Chapter III discusses three MM FORTRAN programs used 
to determine the current on the loop. The entire domain solution (due to R. F. 
Harrington), which uses complex Fourier modes (HARLOOP) is considered to be the 
most accurate and therefore serx'es as a baseline for evaluating the other solutions. A 
second program that uses linear segments (LOOPSCAT) is compared to a third program 
that uses curved segments (CURVENEW). Chapter IV discusses the results obtained by 
the three methods and presents some guidelines in choosing an optimum solution method 
for a given antenna or scattering problem. 
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n. THE THIN WIRE INTEGRAL EQUATION 



A. DERIVATION OF THE THIN- WIRE ELECTRIC FIELD INTEGRAL 

EQUATION 

In this chapter, the integral equation for the current on a thin wire will be 
developed. Time-harmonic field quantities are assumed throughout. Phasor quantities are 
used with the d"' dependency suppressed. 

Referring to the thin wire geometry of Figure 2, the origin is point 0, the location 
of a source point is given by the vector r' and an observation point by the vector r. The 
wire radius, a, is considered constant over the length of the wire. The vector I is 
everywhere parallel to the surface of the wire. Since the sum of the tangential 
components of the incident and scattered electric field must vanish at the surface of a 
perfect electric conductor, the boundary condition is, 

« X (£'+£0 = 0 ( 2 . 1 ) 

If the radius of the wire is small compared to the wavelength of the excitation, the 
surface current density, J,, can be considered constant around the circumference of the 
wire and directed along its axis. The excitation field can be either an incident wave or 
an impressed voltage. The scattered field is the field due to the current on the conductor 
induced by the excitation field. 

The wave equation in terms of the vector potential A is given by 
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V ^A + pU = 



( 2 . 2 ) 



where ^=2-kI\. The solution to equation (2.2) is 



A = g(r,r')dS' 

s 



(2.3) 



471'^ |r-r' 



where the integration is over the primed (source) coordinates. The Green’s function, 
g(r,r') is defined as, 



g(ry) = 



e-mr-r'\ 

4Tz\r-r'\ 



(2.4) 



The expression for the scattered electric field in terms of A is. 
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£* = -ywA V(V-A) 

(i)^e 



(2.5) 



Applying the boundary condition of equation (2.1), 



on S 



( 2 . 6 ) 



o)pe 



Substitution of A in equation (2.3) into equation (2.6) gives, 



^'un = Jo)nfj,g(r,r')dS'+-^V V 
s' 






on S. (2.7) 



Call the second term on the right side of equation (2.7) V, and assume the medium to 
be homogeneous, 



V = V 



V • ^xfj^g(ry)dS' 

s' 



= V 



‘ (-^, g(r,r'))dS' 
s' 



(2.8) 



The vector identity for the divergence of a scalar u times a vector v is, 



V-(wv) = Vwv+w(V-v) 



(2.9) 



Applying this identity to equation (2.8) gives 



V = V 



\^fc7g(.ry))'jys' . 



( 2 . 10 ) 



It can be shown that Vg(r,r') = -Vg(r,r') [Ref. 2]. Using this in equation (2.10) 
and applying the identity of equation (2.9) again yields, 
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y = -V 






( 2 . 11 ) 



=V n/g(r,r')(V'-J,)d5'-^/v'. {g{r,r')J ,)dS 



where V is the del operator defined with respect to the primed coordinates. The second 



integral on the right side of equation (2.11) is equal to zero by the surface divergence 
theorem [Ref. 3]. Thus, V simplifies to 



Substitution of equation (2.12) into equation (2.7) gives an integral equation for 



Equation (2.14) is a form of the Electric Field Integral Equation (EFIE). The unknown 
quantity to be solved for is J,. 

B. SOLUTION OF THE EHE USING MM 

The method of moments (MM) technique can be used to solve for J, by expanding 
it into a series of basis functions, Jj, 




( 2 . 12 ) 




(2.13) 



which may be expressed more compactly as 




(2.14) 
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(2.15) 



t 

where the C, are complex constants to be determined. Substitution of equation (2.15) into 
2.14 gives 



= EC,/[yuny,g(r,r')*-L(V'.J,)Vs(r.r') 

w i. 



dS' 



(2.16) 



To generate the required N equations to solve for the N unknowns, define a suitable 
weighting function \\\, and take the inner product of with both sides of equation 
2.16. The inner product is defined such that it satisfies 



<H’,V> = <V, H’> 

<<x/+yv,w> = a</,H’> +y<v,w> 
<v*,v> >0 if V * 0 
<v’,v> = 0 if V = 0 

[Ref 4], Choose the following inner product: 

s 

where * signifies the complex conjugate. This leads to 



(2.17) 



(2.18) 



fw,-E‘^dS = fji^nW.-ZcJ 



S, 



G)e 



dS'dS . (2.19) 



Interchanging the order of summation and integration and applying the surface divergence 
theorem yields for the right hand side of equation (2.19) 
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J^cj f\joi^i(W^-J,)g(r,r')--^(V'-J,)(V- W^)g(r,r') dS'dS . (2.20) 

• = 1 5 5 . 



h! 



By making the following definitions, 







( 2 . 21 ) 



2.20 and 2.21 can be written in matrix form, 



[K] = [Z][C] 



( 2 . 22 ) 



where [V], [Z] and [C] are called the generalized voltage, impedance and current 
matrices, respectively. The unknown [C] may be solved by an appropriate matrix 
inversion algorithm. Symbolically, 



The generalized current matrix elements are the weighting coefficients in the summation 
of equation (2.15). The current J, is computed from equation (2.15) and the scattered 
field is calculated using this current in the radiation integrals. It should be noted that [V], 
[Z], and [C] have units of volts, ohms, and amperes, but are not unique. In general, they 
depend on the choice of basis and weighting functions. However, the current will 
converge to the same numerical value as the number of basis functions are increased, 
provided the solutions are implemented correctly. 



[C] = [Z]-’[F] . 



(2.23) 
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C. SPECIALIZATION OF THE EFIE TO A CIRCULAR LOOP USING 



CONTORMAL SUBSECTIONS. 

The MM procedure will now be applied to a circular loop in the X-Y plane as 
illustrated in Figure 3. The loop is an ideal test geometry to study the characteristics of 
a MM solution using conformal subsections. It is a relatively simple geometry and other 
accurate solution methods are available to evaluate the performance of conformal 
subsections. The loop has a radius Tq and is divided into N conformal segments. In this 
case a conformal segment is a circular arc. The arclength of the i'*‘ segment is, 



The basis functions, J;, of equation (2.15) are chosen to be overlapping triangles. 




(2.24) 



2 




X 



Figure 3. Thin Wire Loop Geometry 
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(2.25) 



0; elsewhere 



Triangular basis functions are chosen because they are a more accurate representation of 
the current than a pulse basis function since the current is continuous everywhere along 
the wire, and they are relatively easy to deal with analytically. Balanis [Ref. 5] 
states that increasing the basis function complexity beyond triangles may not be 
warranted by the additional improvement in convergence rate. The triangular basis 
functions span two segments and overlap as shown in Figure 3. Therefore, the resultant 
current will be piecewise linear. The weighting (or testing) functions in equation 
(2.19) are chosen such that W^=Ji; (Galerkin’s method). Wang [Ref. 6] states 
that Galerkin’s method provides numerical results which are more accurate than other 
testing methods under similar computational constraints. Substitution of the above 
weighting and basis functions in the expression for Zj, in equation (2.21) yields 




(2.26) 



Using the following relations. 
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/•/' = 4 )- 4 )' = cos(4>-4)') 
/ = r„4); I' = r,4)' 



= iTzar^d^', dS' = 27tarQ<f({)'' 
r,(4)) = r,(ro<t)) = T,{1) 

r (//) = lr^.(4)) 



n = 



N 



— ; P = 

e 



(2.27) 



equation (2.26) may be written as 



Z.V = 



2 n 

'•ovPn 



4ti 



// 



7'i(4>)7'/4'0cos((t. - 4>0 - ^r/(<J)07-/(<|)) 



f‘rl 



(2.28) 



l-RI 






From Figure 4 and the law of cosines, | R | is given by, 



|Rp = IrM -cos(4)-4.')]-fl' 



(2.29) 




Figure 4. Geometry for Determining R. 
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and using the trigonometric identity, 




(2.30) 



results in, 




(2.31) 



By choosing the test (unprimed) points at the center of the wire and the source points on 
the surface of the wire, the singularities along the diagonal of [Z] ai <t>=4)' w'here r=r' 
are avoided. The technique used to calculate | R | is discussed further in Chapters III 
and IV. 

The voltage elements, given in equation (2.21) become 

= '-o f . (2-32) 

'•o't’t 

The incident field, E', for the purpose of this study, will be a plane wave. Figure 5 
shows the direction of incidence of the plane w'ave in spherical polar angles 6 = Q and 
4>=^ measured from the Z and X axes, respectively. E‘ can be 6 or 4> polarized. 
Referring to Figure 5, for 0 polarization, the component of E tangential to the loop, is 



= EqCos© sin(<I>-4)) . 



(2.33) 



Similarly, for 4> polarization. 
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£^6 = £^cos(<I)-4>) . 

The component of the phase vector, /3, parallel to r is 



P - r = sin0 cos(^>-4>) . 



(2.34) 



(2.35) 



Equations (2.30) and (2.32) combine with 2.29 to give 






Vek = r^Eg cos9 J 7^(4>)sin(0-4>)e 



-;prQsin0 cos(<P-i>) 



d<t> . 



(2.36) 
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A similar expression for a <i> directed incident field is 



. ( 2 . 37 ) 

■t. 

The computer coding of the solution for the thin wire loop using the equations 
developed above is described in the next chapter. 
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m. COMPUTER CODES FOR THE THIN WIRE LOOP 



In this section, the FORTRAN program for a thin wire loop using curved segments 
is discussed. The results are presented and compared to similar solutions using straight 
subsections and Fourier modes. 

A. DESCRIPTION OF THE CODES 

Table 1 summarizes the organization of the three programs. Computer listings are 
given in Appendix B and equations from Chapter II will be referenced with a "2." 
preceding the equation number. The FORTRAN source code for the conformal 
subsections is named CURVENEW, and the codes for the straight subsections and the 
Fourier mode solution are named LOOPSCAT and HARLOOP respectively. 

CURVENEW, LOOPSCAT and HARLOOP are functionally similar. Each 
calculates the loop geometry based on the segment size, loop radius, and wire radius, and 
each uses Gaussian quadrature for numerical integration. CURVENEW computes the 
impedance matrix, [Z], in subroutine ZCURVED from the loop geometry of Figure 5 
using equation (2.28). LOOPSCAT uses a similar formulation applied to straight 
segments in subroutine ZMATWW'. HARLOOP computes [Z] using the equations 
developed in reference [7]. CURVENEW computes the excitation vector, [V], using 
equation (2.36) or (2.37) in subroutine CURVEW. LOOPSCAT uses a similar 
formulation for straight subsections in subroutine PLANEW. HARLOOP computes [V] 
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using the equations in reference [7] in subroutine PLANEW. In CURVENEW, all 
integrals are evaluated numerically and symmetry of the impedance elements is used to 
fill the [Z] matrix and reduce the number of numerical calculations. Two formulations 
were investigated with LOOPSCAT: One using a delta function approximation to 
evaluate one of the double integrals in equation (2.21) and the other using Gaussian 
quadrature for both integrations. CURVENEW does all numerical integrations using 
Gaussian quadrature. Matrix symmetry is also used in LOOPSCAT to reduce the number 



TABLE 1. FUNCTIONAL SUMMARY OF PROGRAMS 



Program- > 


CURVENEW 


LOOPSCAT 


HARLOOP 




(Curx’ed 


(Straight 


(Fourier 


Function 

1 


Segments) 


Segments) 


modes) 


Read Input 


Lines 25-38 


Lines 15-27 


Lines 19-30 


Parameters 








Establish 


Lines 45-71 


Lines 53-76 


Lines 44-54 


Loop 








Geometry 








Calculate 


Subroutine 


Subroutine 


Subroutine 


[Z] 


ZCURVED 


ZMATWVV 


ZMATWU 


Calculate 


Subroutine 


Subroutine 


Subroutine 


[V] 


CURVEW 


PLANEW 


PLANEW 


Solve 


Subroutines 


Subroutines 


Subroutines 


System 


DECOMP and 


DECOMP and 


DECOMP and 


[C] = [Z]'[V] 


SOLVE 


SOLVE 


SOLVE 1 

! 

1 


Calculate E* 


Lines 140-180 


Lines 160-197 


Lines 93-137 
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of calculations for [Z], Computation of [C] is performed in subroutines DECOMP and 
SOLVE [Ref. 3], which solve the system of equations using Gaussian elimination. 
Subroutines DECOMP and SOLVE are common to all three programs. The subroutines 
ZCURVED and CURVEW are discussed in more detail in the next section. 

1. Loop Geometry 

To generate the loop geometry an initial estimate of the desired circular arc 
length, Al, is provided. This is used to estimate an angular increment, A<f), 

A<}) = — . (3.1) 

From this estimate, the number of generating points is calculated by. 



N = Int 




and A4> is recalculated using. 



(3.2) 



A(J) = — (3.3) 

N 

to ensure that A({) is such that N segments fill exactly 2x radians. The loop points P, and 
P2 are coincident with P^.j and Pn_2 so that the current is continuous around the loop. 

2. Subroutines ZCURVED and ZMATVVW 

Subroutines ZCURVED and ZMATWW take advantage of the symmetry that 
exists on the loop. For all basis functions, the self impedance terms are equal. The 
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remainder of the matrix is filled with impedance values that repeat in a cyclic manner. 
Mathematically, 



^11 


- Z 22 


- Z 33 - 


= 7 

— ^NN 


II 


^23 = 


^34 =•• 


= z 

• N 


^13 = 


^24 = 


Z 35 


=r 7 

• ^N-2 N 



7 =7 

■^1 N-l ^2 N 

The elements along any diagonal of [Z] are equal and the lower off-diagonal elements 
are the mirror image of the upper diagonals. Thus [Z] is a symmetrical Toeplitz matrix. 
Therefore, computation of the first row of [Z] provides enough information to fill the 
entire matrix. 

Because of the Green’s function in the integrand for the impedance elements, the 
numerical treatment of the self term is very important. To optimize the convergence rate 
and accuracy of CURVENEW and LOOPSCAT, several different approaches are used 
to evaluate | R | near the singularity point where r = r'. In the first method, the 
observation point is chosen along the axis of the wire and the source point along the 
surface for all i,k, giving | R | =a at (/> = (/)' (equation (2.31)). For the second method, 
both the obser\'ation point and the source point are chosen along the axis of the wire 
except on the segment i=k where the value of 4 > at the midpoint is chosen on the axis 
of the wire, with r' on the surface of the wire. Finally, both the observation point and 
the source point are chosen along the axis of the wire, except on the segment i=k, where 
r is chosen along the axis, and r' is chosen on the surface. Choosing the source point and 
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observation point as in the first case gives the most accurate results, but only slightly 
more accurate than the third case. Case two is accurate for small segment sizes but is 
inaccurate for lau’ger segment sizes. Case three was selected for both CURVENEW and 
LOOPSCAT because it is only slightly less accurate than case one, and required fewer 
lines of code. 

ZCURVED calculates the impedance elements of the first row of [Z] by breaking 
the integral in equation (2.28) into four parts. For example, in the first row of [Z], Z,;, 
is calculated by summing contributions from the following four regions of integration 
(Figure 6): 

1 . A double integration along the positive slope of the T, over segment 1 and the 
positive slope of T, over segment i. 

2. A double integration along the negative slope of the T, over segment 2, and the 
positive slope of Tj over segment i. 

3. A double integration along the positive slope of the T, over segment 1, and the 
negative slope of T, over segment i + 1 . 

4. A double integration along the negative slope of the T, over segment 2, and the 
negative slope of T, over segment i + 1. 

The integrations are computed in a similar manner for the derivatives of the basis 
functions, T,' and Tj', over the same subsections. A similar procedure is used for the 
straight subsections in subroutine ZMATWW to calculate the impedance elements. 

The numerical integrations are performed using Gaussian quadrature, with 
the number of points per subsection specified as an input parameter to program 
GAUSWGT, which computes the Legendre polynomial coefficients for a specified 
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number of integration points and writes them to file OUTGLEG. Gaussian quadrature 
w’as chosen because it requires fewer function evaluations than other methods for a given 
accuracy and does not require equal interval samples [Ref. 8]. CURVENEW, 
LOOPSCAT and HARLOOP read the coefficients from file OUTGLEG. The number 
of integration points per wavelength was varied to optimize convergence of LOOPSCAT 
and CURVENEW and is discussed in more detail in Chapter IV. The excitation vector 
[V] is calculated from equation (2.36) or (2.37) in subroutine CURVEW of 
CURVENEW. 

3. Execution Time 

Analysis of the nested DO loop structure of subroutine ZCURVED of 
program CURVENEW indicates that the total execution time of ZCURVED can be 
represented by, 

K = 'X,N,Nl (3.5) 



Tl 



Figure 6. Loop Geometry for Impedance Integrations 
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where N<. is the number of curved segments (from equation (3.2)) is the number of 
Gaussian integration constants per curved segment and is a constant. Execution time 
of the subroutines DECOMP and SOLVE, common to both CURVENEW and 
LOOPSCAT, can be represented similarly by [Ref. 9], 

where N is the number of segments. The excitation subroutine CURVEW execution time 
and field integrations are given by, 

( 3 - 7 ) 

where again y and are constants. Assuming that the execution time of the rest of the 
program is negligible, the total execution time of CURVENEW is 






A similar expression for LOOPSCAT uses the subscript 1, 



(3.8) 



T, = • (3-9) 

Run times were recorded for various values of N^, N, and Ng using an IBM PC/AT with 
a math coprocessor and the coefficients for CURVENEW are found to be 7=0.000156, 
0^=0.0230, and ^^=0.0222. The coefficients for LOOPSCAT are a|=0.0132 and 
^1=0.0252. For the moment, assume that the number of Gaussian integration points per 
wavelength, Ng, is held constant for both CURVENEW and LOOPSCAT. The number 
of integration points on a segment is. 
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and the number of segments is, 



2nro 

= 5 

A/ 



(3.11) 



Similar expressions may be written for straight subsections. Combining equations (3.8), 
3.10, and 3.11 gives 



^2 



N. 



g 



(3.12) 



The ratio of T<, to T, is given by. 



TJT, = 



An^rlNlaJ N^ + yNl + 2 n 
4-n^rlNlaJ N^ + yNf 



(3.13) 



Equation (3.13) will be used in Chapter IV to compare the execution times of 



CURVENEW and LOOPSCAT. 



CALCULATED DATA FOR THE LOOP 



The convergence of the MM solutions for both the current and electric field for 
circular loops of various dimensions are presented for both linear and circular 
polarizations. The convergence is shown to depend on the segment size and number of 
integration points, as well as excitation conditions (incidence direction and polarization). 
Representative plots are presented within the chapter, and additional plots are given in 
Appendix B. 

A. CONVERGENCE OF HARLOOP 

The Fourier mode solution, HARLOOP, was tested for convergence with respect 
to incidence angle, number of modes, and number of integration constants to establish 
a baseline for comparison to CURVENEW and LOOPSCAT. HARLOOP was chosen as 
a baseline because the sinusoidal basis functions match the physical behavior of the 
current on the loop, and thus the current series converges rapidly. This is illustrated in 
Figure 7 for a 0.5 X radius loop with a plane wave incident at an angle of 40 degrees. 

Oscillation of the current as a function of <i> becomes more rapid as 0 is increased, 
because the phase of the incident field over the loop varies as sin 0 (equation (2.35)). 
For a 6 polarized linear wave incident in the 4>=0 plane, the current is always zero at 
4>=0 and 180 degrees, where E’ is cross-polarized with respect to the axis of the wire. 
For 6 polarized incident waves, maxima occur at <t>=90 and 270 degrees, where E is 
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X 1 0-3 



Current Mognitude v« Loop angle 




Figure 7. Convergence of the Current in the Complex Exponential Solution 

parallel to the axis of the wire. The overall current amplitude decreases for © 
approaching 90 degrees, as expected, since the loop’s projected area is small as viewed 
by the incident wave. For 4 > polarized incident waves, the minima occur at 4>=90 and 
270 degrees and maxima at 4>=0 and 180 degrees. The currents do not vanish for 0 
approaching 90 degrees because the loop is parallel to the 4 > polarized incident field. 
Circularly polarized incident waves give a constant magnitude current at normal 
incidence, and oscillations increase with 0. For 0=90 degrees, the circular and </> 
polarization responses are identical. 

HARLOOP is also found to be in agreement with measurements taken on the echo 
area of wire loops at normal incidence [Ref. 10]. The plot of Figure 8 gives the 
echo area (a/X^) versus ro for varying wire radius using HARLOOP. Measured data is 
indicated by the ’ + ’ sign. 
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B. CONVERGENCE OF THE CURRENT EXPANSION 



Having established HARLOOP as a baseline for comparison, convergence of the 
curved subsection program CURVENEW and the linear subsection program LOOPSCAT 
was evaluated. The plots of Figures 9 through 14 give the current on the loop as a 
function of loop angle, 4 >, for varying 0, loop radius, and incident wave polarization. 
The number of integration points per wavelength, Ng, is held constant at 320 in 
LOOPSCAT and CURVENEW. This number was chosen empirically to give a 
converged current within five to ten percent RMS. The RMS error is defined relative to 
the Fourier mode solution. The HARLOOP current is plotted with the solid line, those 
of CURVENEW are plotted with the " + " sign, and those of LOOPSCAT with the "o". 
The wire radius is 0.001 X for these calculations. 
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Figure 8. Backscatter Echo Area for a Loop with varying Radius at Normal Incidence 
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Figure 10. Phase of the Current on a 0.1 X Radius Loop, Normal 
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Figure 11. Magnitude of the Current on a 0.1 X Radius Loop, Incidence 
Angle=40 deg, Circular Polarization (+ =curved; o =linear) 



Current Phose v* Loop ongle 




Figure 12. Phase of the Current on a 0.1 X Radius Loop, Incidence 
Angle=40 deg., Circular Polarization (+ =curved; o =linear) 
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Figure 13. Magnitude of the Current on a 0.1 X Radius Loop, Incidence 
Angle=40 deg., Theta Polarization (+ =curved; o =linear) 
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Figure 14. Phase of the Current on a 0.5 X Radius Loop, Incidence 
Angle=40 deg., Theta Polarization (+ =curved; o =linear) 
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Representative plots of the normalized root mean squared error are given in Figures 15 
through 22. For set values of N and Ng, CURVENEW converges faster than 
LOOPSCAT in most cases. The error difference is most pronounced for loop 
circumferences on the order of a wavelength or less. From the plots, for ro=0. 1 X, 
CURVENEW converges to less than 10 percent error for segment sizes ranging from 
0.02 to 0.2 wavelengths (N(.=32 to N<,=3). LOOPSCAT converges to within 10 percent 
error for segment sizes less than approximately 0.06 X (N, > 10) but gives errors of 30 
to 40 percent for a segment size of 0.2 wavelengths. There is no improvement using 
curved subsections on larger loops for off-axis incidence waves, but the curved 
subsections give small improvements for large loops at normal incidence. This is 
expected in view of the behavior of the current on the loop. For linear polarization the 
current abruptly flips polarity from one side of the loop to the other. For circular 
polarization, the amplitude is constant, but the phase is linear. Both of these conditions 
can be represented accurately by a few triangles if the impedance and excitation integrals 
are evaluated precisely on the loop contour (see Figure 23). 

Plots of the magnitude of the backscattered E' versus incidence angle for varying 
segment size, Ng, and ro are given in Figures 24 and 25. As with the currents, the 
backscattered field converges more rapidly for CURVENEW than LOOPSCAT in most 
cases, with the greatest difference for small radius loops and angles near the maximum 
values of | E* | . 
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Figure 15. Error in the Current Magnitude for a 0.1 X Radius Loop, 
Normal Incidence, Circular Polarization (+ =curved; o =linear) 
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Figure 16. Error in the Current Magnitude for a 0.1 X Radius Loop, 
Incidence Angle=40 deg.. Circular Polarization (+ =cur\'ed; o =linear) 
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Figure 17. Error in the Current Magnitude for a 0.5 X Radius Loop, 
Normal Incidence, Circular Polarization (+ =cur\ed; o =linear) 
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Figure 18. Error in the Current Magnitude for a 0.5 X Radius Loop, 
Incidence Angle=40 deg.. Circular Polarization (+ =curved; o =linear) 
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Figure 19. Error in the Current Magnitude for a 0.1 X Radius Loop, 
Normal Incidence, Theta Polarization (+ =curved; o =linear) 
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Figure 20. Error in the Current Magnitude for a 0.1 X Radius Loop, 
Angle of Incidence =40 deg., Theta Polarization (+ =curved; o =linear) 
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Figure 21. Error in the Current Magnitude for a 0.5 X Radius Loop, 
Normal Incidence, Theta Polarization (+ =curved; o =linear) 
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Figure 22. Error in the Current Magnitude for a 0.5 X Radius Loop, 
Incidence Angle=40 deg., Theta Polarization (+ =curved; o =linear) 
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Figure 23. A Representation of a Sinusoidal Current with Two Basis Functions (top) and 
a Constant Current as a Superposition of Triangles (bottom) 
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Convergence of Bockscottered Field 




Figure 24. Backscattered Electric Field Intensity for varying Angles of 
Incidence, 0.1 X Radius Loop, Theta Polarization (+ =curved; o =linear) 
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Figure 25. Backscattered Electric Field Intensity for varying Angles of 
Incidence, 0.5 X Radius Loop, Theta Polarization (+ =curved; o =linear) 
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C. COMPARISON OF EXECUTION TIME AND MEMORY REQUIREMENTS 



The average run time for a given loop radius and convergence error is greater for 
CURVENEW than for LOOPSCAT due to the dependence in equation (3.13) and 
relative magnitudes of the coefficients o; and y. The plot of equation (3.13) in Figure 26 
illustrates the ratio of run times of CURVENEW and HARLOOP versus N, for N(.=4, 




0 20 40 60 80 100 120 140 160 180 200 

Number of Lineor Subsections 

Figure 26. Ratio of Execution Times for CURVENEW and LOOPSCAT for Fixed 
Curved Segment Lengths and with Varying Linear Segment Lengths 

8, 32 and Ng=320 for a 0.1 X loop. The run time of CURVENEW is less than 

HARLOOP for T^/T, < 1. From the plot, the "break even" points are approximately 

N,= 115, 90, 52 for Ne=4, 8, 32. For N, less than these values, the integration 

subroutine ZCURVED is the determining factor in the run time; for N, greater than these 

values, Gausssian elimination subroutines DECOMP and SOLVE are the determining 

factors. The increased number of integrations per segment completely offset the time 
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savings of a reduced [Z] matrix in CURVENEW. As mentioned in Chapter III, a delta 
function approximation for the outer integration of the impedance integral of equation 
(2.28) was investigated. This reduces the exponent of Ng to one in equation (3.13), but 
many more segments are required for a given accuracy. A more efficient integration 
scheme using a large number of integration points in the vicinity of <t>=<j>' and fewer 
integrations elsewhere may reduce the execution time. 

The savings in computer memory is significant for CURVENEW, since the number 
of matrix elements is on the order of N\ For a given accuracy for a 0. 1 X loop, the ratio 
of the number of elements required for curved and linear subsections, (Nj/Ni)^ is on the 
order of 0.1 to 0.2. This is a reduction of 80 to 90 percent. Although the memory 
requirements for the small loops considered here are not prohibitive for piecewise linear 
segments, greater memory savings will be realized for larger geometries comprised of 
many small features. 
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V. CONCLUSIONS 



The use of conformal subdomain basis functions (curved subsections) to represent 
the current on a thin curved wire was investigated by solving the thin wire electric field 
integral equation using the method of moments. A solution using triangular basis 
functions was computer coded in FORTRAN and validated by comparing it to measured 
data and the results of two other method of moments solutions (LOOPSCAT and 
HARLOOP). The effect of varying loop radius, segment size, number of integration 
points and incident wave parameters on the accuracy and rate of convergence of the 
current expansion and backscattered field was investigated. 

For small loops with circumferences on the order of a wavelength, the number of 
segments required to converge to a given accuracy with the curved segments was as 
small as 20 percent of the number of linear segments required to converge to the same 
accuracy (see Table 2). From computed data, it was determined that the greatest 
reduction in the number of unknowns for curved subsections occurs for geometries where 
the current amplitude variations over the surface are small and the phase variations are 
small or linear. As mentioned in Chapter I, the general rule of thumb for the length of 
one segment is 0.05 X to 0.1 X, which corresponds to a phase variation of 20 to 40 
degrees. This restriction in phase variation is the driving factor when choosing the 
segment size for geometrical features having radii of curvature of approximately X/2 or 
larger. A piecewise linear segment is small enough to represent the geometry accurately 
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in this case. The phase restriction applies to curved segments as well, but the curved 
segments conform exactly to the wire, and hence for small loops there is no sacrifice in 
geometrical accuracy by choosing segment sizes of approximately 0.05 X or 20 electrical 
degrees. As the loop becomes larger, or the wavelength becomes smaller, the curved and 
straight subsection solutions become equivalent. The greatest advantage in using curved 
subsections to reduce the number of segments is for electrically small structures where 
small linear segments are required simply to reproduce the wire shape. 

Although the number of segments was greatly reduced using conformal subsections, 
the execution time was increased due to the increased number of integration points p>er 
segment required for acceptable accuracy. To reduce the integration time, it is suggested 
that the number of integration points per wavelength be varied from a large number when 
evaluating the self term, to fewer points away from the self term. For certain geometries, 
symmetry could also be used to reduce the integration time. 

To avoid singularities, the MM testing procedure was performed along the axis of 
the wire and the current constrained to the surface of the wire. A delta function 
approximation in the impedance integrations was found to reduce the time required to 
compute the impedance matrix, but as expected, required more segments for a given 
convergence accuracy. 

A disadvantage in the formulation of CURVENEW is that an analytic expression 
for the curve is needed to perform the integrations for the impedance matrix and the 
excitation vector. The expressions will change each time a new curve is analyzed, and 
consequently, considerable effort will be required to modify the code every time a change 
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is made in the geometry. Programs that use linear segments are more flexible because 
they require only the coordinates of the points along the wire axis to generate the 
integration points. 

The next logical step in testing the effectiveness of conformal subdomain basis 
functions is to formulate the solution for an equiangular spiral wire. The equiangular 
spiral is used in broadband antennas and has a simple mathematical form. For 
geometrical accuracy, the segment size in the piecewise linear formulation will be much 
smaller than 20 electrical degrees near the center of the spiral. Equal length conformal 



TABLE 2. COMPARISON OF CURVENEW AND LOOPSCAT 





Curved 

Subsections 


Linear 

Subsections 


Loop Radius 


O.IX 


0.5X 


O.IX 


0.5a 


Number of Segments for 10% 
RMS Error. Normal Incidence 


3 


5 


21 


5 


Number of Segments for 10% 
RMS Error, Off Axis Incidence 


3 


26 


10 


26 


Execution Time', 10% RMS 
Error, Normal Incidence 


314 s 


4669 s 


32 s 


2691 s 


Execution Time', 10% RMS 
Error, Off Axis Incidence 


314 s 


918 s 


59 s 


540 s 


[Z] Matrix Size, 10 % RMS 
Error, Normal Incidence 


1 9 


25 


441 


25 


[Z] Matrix Size, 10 % RMS 
Error, Off Axis Incidence 


9 


676 


100 


676 1 



* Execution time measured with an IBM PC/AT. 
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segments may be used along the spiral arms, and it is anticipated that the number of 
required segments will be substantially reduced. 
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APPENDIX A 



COMPUTER CODES 



C MAIN PROGRAM *CURVNEW.FOR* 

C MORE NUMERICALLY EFFICIENT THAN CURVSUB.F •*** 

C PLANE WAVE SCATTERING FROM A CIRCULAR LOOP IN THE Z PLANE. 

C METHOD OF MOMENTS WTTH CURVED SUBSECTIONS 
C **** RADAR CROSS SECTION CALCULATION ***** 

C 

COMPLEX Z(25000),B(500),C(500),BP(500),CP(500),EX,EC 

COMPLEX ET,EP,UC,RW(500),U0 

DIMENSION ECP(400),IPS(500),ANG(400),EDP(400) 

DIMENSION XG(128),AG(128),PH(500),DEL(500),PA(500) 

DIMENSION ECV(5Cf0),EXV(500),PHC(500),PHX(500) 

DATA PI/3.14159265/ 

DATA IPRINT/1/,ITEST/1/ 

RAD = PI/180. 

ECX = 0. 

BK = 2.*PI 
ETA = 377. 

U0 = (0.,0.) 

UC = (0.,-l.)*ETA*BK/4./PI 
CONST=16.*PI*“3 
PHIRD = 0. 

READ INPUT AND PROGRAM CONTROL PARAMETERS 
OPEN( 1 .FILE = ‘PAR AM LST. DAT') 

READ(1,*)ANGLE,DT, START, STOP, AW, RO,SEG,HARR,GCONST,XLOC,XIPOL 
LOC = INT(XLOC) 

IPOL = INT(XIPOL) 

CLOSE(l) 

READ GAUSSIAN CONSTANTS 

OPEN(2,FILE= ‘OUTGLEG’) 

READ(2,*) NG 
DO 1 1=1, NG 
READ(2,*) XG(I),AG(I) 

1 CONTINUE 
C 

C OUTCURV IS THE FILE THAT THE SCATTERED FIELD DATA IS WRITTEN TO 
C ICURVOUT IS THE FILE THAT THE LOOP CURRENT DATA IS WRITTEN TO 
C 
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43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 



0PEN(8,F1LE = ’OURCURV.DAT’) 

0PEN(7,F1LE= ’ICURVOUT.DAT') 

C 

C GENERATE THE LOOP POINTS 

C 

C 

C CHOOSE THE NUMBER OF POINTS BASED ON THE VALUE OF SEG 
C 

C AW = AW*R0 

DPHI = SEG/R0 
NP=1NT(2.*PI/DPHI) + 1 
DPHI = 360. /FLOAT(NP) 

PH(1) = 0. 

DO 10 1 = 2,NP+1 
PH(I) = FLOAT(I- 1 )*DPHI*RAD 
DEL(I-1) = (PH(I)-PH(I-1)) 

P A(I- 1 ) = (PH(I) + PH(I- 1 ))/2 . 

10 CONTINUE 
NP=NP+2 

C 

C OVERLAP THE ENDS SO THAT CURRENT WILL BE CONTINUOUS ON THE LOOP 
C 

PH(NP) = BK + PH(2) 

DEL(NP-1) = DEL(1) 

PA(NP-1)=BK + PA(1) 

MT=NP-2 
DO 52 1=1, NP 
XHB = R0*COS(PH(D) 

YHB = R0*SIN(PH(I)) 

52 CONTINUE 

WRITE(6,*) ’GEOMETRY DEFINED’ 

IF(ITEST.EQ.O) GO TO 98 
C 

C DEFINE DIMENSIONS OF THE IMPEDANCE MATRIX BLOCKS 
C 

WRITE(6,*) ’NP,MT=’,NP,MT 
C 

C COMPUTE IMPEDANCE MATRIX ELEMENTS 
C 

CALL ZCURVED(NP,RO,PH,DEL,PA,NG,XG,AG,AW,Z,ZMAX) 

DO II I=1,MT 
CZ=CABS(Z(I)) 

AZ = ATAN2(AIMAG(Z(I)),REAL(Z(I)) + I.E-8)/RAD 

11 CONTINUE 

WR1TE(6,*) ’WIRE IMPEDANCE COMPUTED’ 

C 

C PERFORM LU DECOMPOSITION 
C 

CALL DECOMP(MT,IPS,Z) 

WRITE(6,*) ’Z DECOMPOSED’ 
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93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

no 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 
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127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 



C 

C BEGIN FIELD CALCULATIONS 
C 

98 PHRO=PHIRD*RAD 
IT=INT((STOP-START)/DT)+ 1 
WRITE(7,*) IT,MT,0,0 

DO 500 1 = 1, IT 

THETA = FLOAT(M )*DT + START 

THR=THETA*RAD 

PHR=PHR0 

IF(THETA.LT.180.) GO TO 99 
THR = (360.-THETA)’^RAD 
PHR = PHRO + PI 

99 CONTINUE 
ET=U0 

EP = U0 
C 

C COMPUTE THE EXCITATION VECTOR 
C 

CALL CURVEW(NP,RO,PH,DEL,PA,NG,XG,AG,THR,PHR,RW) 
IF (LOC .EQ. 0) THEN 
C 

C CIRCULAR POLARIZATION IF LOC=l ELSE LINEAR 
C 

IF(IPOL.EQ.l)THEN 

C 

C THETA POLARIZED INCIDENT WAVE (IPOL=l) 

C 

DO 101 L=1,MT 
101 B(L) = RW(L) 

ELSE 

C 

C PHI POLARIZED INCIDENT WAVE (IPOL = 2) 

C 

ENDIF 

IF(ITEST.EQ.O) GO TO 9998 
C 

C PERFORM GAUSSIAN ELIMINATION TO DETERMINE [C] 

C 

CALL SOLVE(MT,IPS,Z,B,C) 

DO 210 L=1,MT 

WRITE(7,*) L,L’^DPHI,CABS(C(L)),ATAN2(REAL(C(L)), 

AIMAG(C(L)))/RAD 

ET=ET + RW(L)*C(L) 

210 EP = EP + RW(L + MT)*C(L) 

ELSE 

C 

C THETA POLARIZED INCIDENT WAVE 
C 

DO 221 L=LMT 
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221 B(L)=RW(L) 

C 

C PHI POLARIZED INCIDENT WAVE 
C PHASE SHIFT FOR CP IS PI/2. 

C 

DO 222 L=1,MT 

222 BP(L)=RW(L+MT)*CEXP(CMPLX(0.,PI/2.)) 

CALL SOLVE(MT,IPS,Z,B,C) 

CALL SOLVE(MT,IPS,Z,BP,CP) 

DO 220 L=1,MT 

WRITE(7 , ♦) L,L*DPHI , C ABS(C(L) + CP(L)) , ATAN2(RE AL( C(L) + CP(L)) , 
♦ AIMAG(C(L) + CP(L)))/RAD 

ET = ET + R W(L) *C(L) + R W(L) *CP(L) 

220 EP= EP + RW(L + MT)*C(L) + RW(L + MT)*CP(L) 

ENDIF 
EC = ET*UC 
EX = EP*UC 
ANG(I) = THETA 
ECV(I) = CABS(EC) 

EXV(I) = CABS(EX) 

ECR=REAL(EC) 

ECI = AIMAG(EC) 

EXR=REAL(EX) 

EXI = A1MAG(EX) 

PHC(I) = ATAN2(EC1,ECR+ 1.E-20)/RAD 
PHX(I) = ATAN2(EXI,EXR+l.E-20)/RAD 
ECX = AMAX1(ECX,ECV(I),EXV(I)) 

500 CONTINUE 

WRITE(6,*) ’EMAX = ’,ECX 
DO 600 I = 1,IT 

ECV(I)=AMAX1(ECV(I),1.E-10) 

EX V(I) = AM AX 1(EXV(I), 1 .E- 1 0) 

ECP(I) = (EC V(I)/ECX)**2 
EDP(I) = (EX V(I)/ECX)**2 
ECP(I) = AM AX1(ECP(I),. 00001) 

EDP(I) = AMAX1(EDP(1), .00001) 

ECP(1)= 10.*ALOG10(ECP(I)) 

EDP(1)= 10.*ALOG10(EDP(I)) 

600 CONTINUE 

SIGM A = (ECX**2)*CABS(UC)/(2.*BK) 

SIGDB=10.*ALOG10(SIGMA) 

WRITE(6,*) ’BACKSCATTER CROSS-SECTION, IN DB = ’,SIGMA,SIGDB 
208 F0RMAT(/,5X,’SIGMAAVAVLSQ = ’,E15.4, 

*/,5X,’ INDB=’,F8.4) 

DO 9000 L= LIT 
WRITE(8,*) ANG(L),ECV(L) 

9000 CONTINUE 
9998 STOP 
END 

SUBROUTINE SOLVE(N,IPS,UL.B,X) 
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SUBROUTINE TO SOLVE SYSTEM OF EQUATIONS WITH COMPLEX 
OEFFICIENTS. 

CALL ’DECOMP' FIRST. (FROM MAUTZ AND HARRINGTON) 

COMPLEX UL(5(XX)0),B(200),X(200),SUM 
DIMENSION IPS(5(X)) 

NP1=N+1 

1P=IPS(1) 

X(1) = B(IP) 

DO 2 1=2, N 
IP = IPS(I) 

IPB = IP 
IM1=I-1 
SUM = (0.,0.) 

DO 1 J=1,IM1 

SUM = SUM + UL(IP)*X(J) 

1 IP = IP + N 

2 X(I) = B(IPB)-SUM 
K2 = N*(N-1) 

IP = IPS(N) + K2 
X(N) = X(N)/UL(IP) 

DO 4 IBACK = 2,N 

I = NP1-IBACK 

K2=K2-N 

IPI = IPS(I) + K2 

IP1=I+1 

SUM = (0.,0.) 

IP = IPI 

DO 3 J = 1P1,N 
IP = IP + N 

3 SUM = SUM + UL(IP)*X(J) 

4 X(I) = (X(I)-SUM)/UL(IPI) 

RETURN 

END 

SUBROUTINE DECOMP(N,IPS,UL) 



C SUBROUTINE TO DECOMPOSE SYSTEM OF EQUATIONS. 
C FROM MAUTZ AND HARRINGTON. 

C 

COMPLEX UL(50000),PIVOT,EM 
DIMENSION SCL(200),IPS(200) 

DO 5 1 = 1, N 
IPS(I) = I 
RN=0. 

J1 = I 

DO 2 J=1,N 

ULM = ABS(RE AL(UL(J 1 ))) + ABS(AJM AG(UL(J 1 ))) 
J1=J1+N 
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IF(RN-ULM) 1,2,2 
1 RN=ULM 
2 CONTINUE 
SCL(1)=1./RN 
5 CONTINUE 
NM1=N-1 
K2 = 0 

DO 17 K=1,NM1 
BIG=0. 

DO 11 I=K,N 
IP=IPS(I) 

IPK=IP + K2 

SIZE = (ABS(REAL(UL(IPK))) + ABS( AIM AG(UL(IPK))))*SCL(IP) 

IF(SIZE-BIG) 11,11,10 

10 BIG = SIZE 
IPV=I 

11 CONTINUE 
IF(IPV-K) 14,15,14 

14 J = IPS(K) 

IPS(K) = IPS(IPV) 

IPS(IPV)=J 

15 KPP = IPS(K) + K2 
PIVOT = UL(KPP) 

KP1=K+1 

DO 16 I = KP1,N 
KP = KPP 
IP = IPS(I) + K2 
EM = -UL(IP)/PIVOT 
18 UL(IP) = -EM 
DO 16 J = KP1,N 
IP = IP + N 
KP=KP+N 

UL(IP)=UL(IP) + EM*UL(KP) 

16 CONTINUE 
K2 = K2 + N 

17 CONTINUE 
RETURN 
END 

SUBROUTINE ZCURVED(NP,RO,PH,DEL,PA,NG,XG,AG,A,Z,ZMAX) 

C 

C IMPEDANCE ELEMENTS FOR CURVED BASIS FUNCTIONS. 

C SPECIFICALLY DERIVED FOR A CIRCULAR LOOP - NOT A GENERAL CASE. 
C 

COMPLEX CEXP,Z(50000),CON,CMPLX,SUMA,SUMB 
COMPLEX U0,SUM1,SUM2,SUM3,SUM4,EXP,ZT(500) 

DIMENSION DEL(500),PH(500),XG(128),AG(128),PA(500) 

C OPEN(2,FILE = ’ZCURV.DAT’) 

ETA = 377. 

ZMAX = 0. 

PI = 3. 14159 
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BK = 2.*PI 
BK2 = BK**2 
U0=(0.,0.) 

CON = (0. . 1 . )*BK*ETA/(4 . *PI)*R0**2 
NT = NP-2 
C 

C COMPUTE Z(1,LQ) = ZT(LQ) 

C 

KQ=1 

Pl=DEL(KQ)/2. 

P2=PA(KQ) 

P3=DEL(KQ + l)/2. 

P4 = PA(KQ+1) 

C 

C DO THE L LOOP 
C 

DO 600 LQ=1,NT 
PPl=DEL(LQ)/2. 

PP2=PA(LQ) 

PP3 = DEL(LQ+l)/2. 

PP4=PA(LQ+1) 

C 

C DO THE PHI INTEGRATION 
C *** FIRST PART FROM PHI(K) TO PHI(K + 1) 

C PHI PRIMED INTEGRATION FOR THE POSITIVE SLOPE OF LQ 
C 

SUMA = U0 
PHA = PA(KQ) 

DO 100 1=1, NG 
PH1 = P1*XG(I) + P2 
TK = (PHI-PH(KQ)VDEL(KQ) 

TKP=1./DEL(KQ)/R0 
SUM1 = U0 
DO 90 J=1,NG 
PHIP = PP1*XG(J) + PP2 
TL=(PHIP-PH(LQ))/DEL(LQ) 

TLP=1./DEL(LQ)/R0 
CC = C0S(PH1-PH1P) 

C 

C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
C 

RR = R0*SQRT(4. *(SIN((PHI-PHIP)/2. ))**2 + ( A/R0)**2) 

EXP=CEXP(CMPLX(0.,-BK*RR))/RR 

SUMl = SUM1 + AG(J)*EXP*(TK*TL*CC-TKP-n'LP/BK2) 

90 CONTINUE 

SUM1 = SUMI*PP1 
C 

C PHI PRIMED INTEGRATION FOR THE NEGATIVE SLOPE OF LQ 
C 
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SUM2 = U0 
DO 80 J = 1,NG 
PH1P = PP3*XG(J) + PP4 
TL= 1 . -(PHIP-PH(LQ + 1 ))/DEL(LQ + 1 ) 

TLP = -1./DEL(LQ+1)/R0 
CC = COS(PHl-PHIP) 

C 

C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
C 

RR=R0*SQRT(4.*(SIN((PHI-PHIP)/2.))*’^ + (A/R0)**2) 

EXP=CEXP(CMPLX(0.,-BK*RR))/RR 

SUM2 = SUM2 + AG(J)+EXP*(TK*TL*CC-TKP*TLP/BK2) 

80 CONTINUE 

SUM2 = SUM2*PP3 

SUMA = SUMA + (SUM1+SUM2)*AG(I) 

100 CONTINUE 

SUMA = SUMA*P1 
C 

c *** SECOND PART FROM PHI(K+ 1) TO PHI(K + 2) 

C PHI PRIMED INTEGRATION FOR THE POSITIVE SLOPE OF LQ 
C 

SUMB=U0 
PHA = PA(KQ+1) 

DO 101 I=1,NG 
PHI = P3*XG(I) + P4 
TK= 1 ,-(PHI-PH(KQ+ 1))/DEL(KQ+ 1) 

TKP = -1./DEL(KQ+1)/R0 
SUM3=U0 
DO 91 J = 1,NG 
PHIP = PP1*XG(J) + PP2 
TL= (PH1P-PH(LQ))/DEL(LQ) 

TLP=1./DEL(LQ)/R0 
CC = COS(PHl-PHIP) 

C 

C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
C 

RR = R0*SQRT(4. *(SIN((PHI-PHIP)/2.))**2 + (A/R0)**2) 

EXP=CEXP(CMPLX(0.,-BK*RR))/RR 

SUM3 = SUM3 + AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 

91 CONTINUE 

SUM3 = SUM3*PP1 
C 

C PHI PRIMED INTEGRATION FOR THE NEGATIVE SLOPE OF LQ 
C 

SUM4 = U0 
DO 81 J = 1,NG 
PHIP = PP3*XG(J) + PP4 
TL= l.-(PHIP-PH(LQ+ 1))/DEL(LQ+ 1) 
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TLP = -1 . /DEL(LQ + 1 )/R0 
CC=COS(PHI-PHIP) 

C 

C COMPUTE THE MAGNITUDE OF R. NOTE THAT R IS COMPUTED FROM THE 
C WIRE AXIS TO THE SURFACE OF THE WIRE 
C 

RR=R0*SQRT(4.*(SIN((PHI-PHIP)/2.))**2+(A/R0)**2) 

EXP = CEXP(CMPLX(0. ,-BK*RR))/RR 

SUM4 = SUM4+AG(J)*EXP*(TK*TL*CC-TKP*TLP/BK2) 

81 CONTINUE 

SUM4 = SUM4*PP3 

SUMB = SUMB +(SUM3 + SUM4)*AG(I) 

101 CONTINUE 

SUMB=SUMB*P3 

ZT(LQ) = CON*(SUM A + SUMB) 

ZM AX = AM AX 1(ZM AX,CABS(ZT(LQ))) 

600 CONTINUE 
ZT(NT) = ZT(2) 

C 

C FILL THE ENTIRE Z MATRIX USING SYMMETRY PROPERTIES 
C [Z] IS A SYMMETRICAL TOEPLITZ MATRIX 
C ROW INDEX, I; COL INDEX, J 
C 

DO 10 1=1, NT 
DO 10 J=1,NT 
K = (1-1)*NT + J 
Z(K) = U0 
U = IABS(I-J) 

IF(IJ.GT.NT) GO TO 10 

IJ1=IJ+1 

Z(K) = ZT(IJ1) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE CURVEW(NP,RO,PH,DEL,PA,NG,XG,AG,THR,PHR,R) 

C 

C PLANE WAVE EXCITATION VECTOR ELEMENTS FOR A LOOP USING CURVED 
C BASIS FUNCTIONS. INCIDENCE DIRECTION IS (THR,PHR). THE WIRE LIES 
C IN THE X-Y PLANE (Z=0). 

C 

COMPLEX U0,R(500),SUM 1 ,SUM2,SUM3,SUM4,CEXP,FF 

COMPLEX GG,CMPLX,SUMT,SUMP 

DIMENSION PH(500),DEL(500),PA(500),XG(128),AG(128) 

MT=NP-2 
U0 = (0.,0.) 

PI = 3. 14159 
BK = 2.*P1 
ST=SIN(THR) 

CT=COS(THR) 

DO 50 IP=1,MT 
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SUM1=U0 
SUM2=U0 
SUM3 = U0 
SUM4=U0 
Pl=DEL(IP)/2. 

P2 = PA(IP) 

P3=DEL(IP+l)/2. 

P4 = PA(1P+1) 

DO 20 1=1, NG 
PHI = P1*XG(I) + P2 
CC = COS(PHR-PHl) 

SS = SIN(PHR-PHI)*CT 

FF = AG(I)*(PHI-PH(IP))/DEL(IP)*CEXP(CMPLX(0.,BK*RO*ST*CC)) 

SUM1=SUM1+CC*FF 

SUM2 = SUM2 + SS*FF 

PH1 = P3*XG(1) + P4 

CC = COS(PHR-PHI) 

SS = SIN(PHR-PHI)*CT 

GG = AG(1)*(1.-(PHI-PH(1P+ 1))/DEL(1P+ 1))*CEXP(CMPLX(0., 

* BK*R0*ST*CC)) 

SUM3 = SUM3 + CC*GG 
SUM4 = SUM4 + SS*GG 
20 CONTINUE 

SUMP= SUM 1 *P1 + SUM3*P3 
SUMT = SUM2*P1 +SUM4*P3 

R-WIRE-THETA IN R(IP) AND R-WIRE-PHI IN R(IP + MT) 

R(IP) = SUMT*R0 
R(IP + MT) = SUMP*R0 
50 CONTINUE 
RETURN 
END 
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C MAIN PROGRAM *LOOP.FOR* 

C PLANE WAVE SCATTERING FROM A CIRCULAR LOOP IN THE Z PLANE. 

C ♦*** RADAR CROSS SECTION CALCULATION 
C 

COMPLEX Z{15000),B{500),C(500),BP(500),CP{500),U 

COMPLEX ET,EP,UC,RW(500),U0 

DIMENSION ECP(400),IPS(500),ANG(400),EDP(400) 

DIMENSION ZH{200),XT(128),AT(128),XH{200),YH(200) 

DIMENSION ECV(400),EXV(400),PHC(400),PHX(400) 

DATA PI/3.14159265/ 

DATA IPRINT/1/ 

READ INPUT AND PROGRAM CONTROL PARAMETERS 
OPEN(l ,FILE = ’PARAMLST.DAT’) 

RE AD(1,*)ANGLE,DT, START, STOP, AW, RB,SEG,HARR,GCONST,XLOC,XIPOL 
LOC = INT(XLOC) 

IPOL = INT(XIPOL) 

CLOSE(l) 

READ GAUSSIAN CONSTANTS 

OPEN(2,FILE= ’OUTGLEG’) 

READ(2,*) NT 
DO 2 I=1,NT 
READ(2,*) XT(I),AT(1) 

2 CONTINUE 
R.AD = PI/180. 

ECX = 0. 

BK = 2.*PI 
ETA = 377. 

U = (0.,1.) 

U0 = (0.,0.) 

UC = -U*ETA*BK/4./PI 
CONST=16.*PI**3 
NT2 = NT/2 
C 

C OUTLOOP IS THE FILE THAT THE SCATTERED FIELD DATA IS WRITTEN TO 
C ISTOUT IS THE FILE THAT THE LOOP CURRENT DATA IS WRITTEN TO 
C 

OPEN(8,FlLE= ’OUTLOOP.DAT’) 

OPEN(7,FILE = TSTOUT.DAT') 

DPH1 = SEG/RB 

NP = INT(2.*PI/DPHI)+1 

DPHI = 360./FLOAT(NP) 

WRITE(6,*) ’DPHI,NP = ’,DPHI,NP 
C 

C GENERATE THE LOOP POINTS. MULTIPLY ALL QUANTITIES BY BK ( = 2*PI) 

C 

C 
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C CHOOSE THE NUMBER OF POINTS BASED ON THE VALUE OF SEG 
C 

AK = AW*BK 

DO 10 I=1,NP+1 

PP = FLOAT(l-l)*DPHl*RAD 

XH(I)=RB*COS(PP)*BK 

YH(I) = RB*S1N(PP)*BK 

ZH(I) = 0. 

10 CONTINUE 
NP=NP+2 
C 

C OVERLAP THE ENDS SO THAT THE CURRENT IS CONTINUOUS 
C 

XH(NP)=XH(2) 

YH(NP) = YH(2) 

ZH(NP) = ZH(2) 

DO 52 I=1,NP 
YHB = YH(I)/BK 
XHB = XH(I)/BK 
DEL = 0. 

IF(I.NE.1)THEN 
DXX = XHB-XH(I-1)/BK 
DYY = YHB-YH(I-1)/BK 
DEL= SQRT(DXX**2 + DYY**2) 

ENDIF 

52 CONTINUE 
C 

C DEFINE DIMENSIONS OF THE IMPEDANCE MATRIX BLOCKS 
C 

MT=NP-2 

WRITE(6,*) ’MT=’,MT 
C 

C COMPUTE IMPEDANCE MATRIX ELEMENTS 
C 

CALL ZMAT\VY\'(1.1,NP,RB,XH,YH,ZH,NT.XT,AT.AK,Z) 

WR1TE(6,*) ’Z COMPUTED' 

IF(IPRINT.EQ.O) THEN 
DO 1010I=1,MT 
CZ = CABS(Z(I)) 

AZ = ATAN2(AlMAG(Z(I)),REAL(Z(I))+l.E-8)/RAD 
WRITE(6,*) 'I, Z= ’,LZ(I), CZ.cz/ZMAX, AZ 
1010 CONTINUE 
ENDIF 
C 

C PERFORM LU DECOMPOSITION 
C 

CALL DECOMP(MT,IPS,Z) 

WRITE(6,*) 'Z DECOMPOSED’ 

C 

C BEGIN FIELD CALCULATIONS. PHI FOR PATTERN CUT (DEGREES) = PHID 
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C 

PHID = 0. 

PHRO = PHID’"RAD 
IT=INT((STOP-START)/DT)-f 1 
WRJTE(7,*) IT,MT,0,0 
DO 500 1=1, IT 

THETA = FLO AT(1 - 1 ) ’*'DT -f START 

THR=THETA*RAD 

PHR=PHR0 

IF(THETA.LT.180.) go to 99 
THR = (360. -THETA) ’►RAD 
PHR = PHR0 + P1 
99 CONTINUE 
ET = U0 
EP=U0 
C 

C COMPUTE THE EXCITATION VECTOR 
C 

CALL PLANEW(NP,XH,YH,ZH,THR,PHR,RW) 

IF (LOC .EQ. 0) THEN 
C 

C CIRCULAR POLARIZATION IF LOC=l ELSE LINEAR 
C 

IF(IPOL.EQ.l)THEN 

C 

C THETA POLARIZED INCIDENT WAVE (IPOL=l) 

C 

DO 101 L=1,MT 

101 B(L) = RW(L) 

ELSE 

C 

C PHI POLARIZED INCIDENT WAVE (IPOL = 2) 

C 

DO 102 L=1,MT 

102 B(L) = RW(L + MT) 

ENDIF 

C 

C PERFORM GAUSSIAN ELIMINATION TO DETERMINE [C] 

C 

CALL SOLVE(MT,IPS,Z,B,C) 

DO 210 L=1,MT 

WRITE(7 , *) L,L*DPHI , CABS(C(L)), ATAN2(REAL(C(L)), 
* AIMAG(C(L)))/RAD 

ET = ET -f (RW(L)/BK)*C(L) 

210 EP=EP-f(RW(L + MT)/BK)*C(L) 

ELSE 

C 

C THETA PLOARIZED INCIDENT WAVE 
C 

DO 221 L=1,MT 
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221 B(L) = RW(L) 

C 

C PHI POLARIZED INCIDENT WAVE 
C PHASE SHIFT FOR CP IS PI/2. 

C 

DO 222 L=1,MT 

222 BP(L) = RW(L+ MT)*CEXP(CMPLX(0. ,PI/2.)) 

CALL SOLVE(MT,IPS,Z,B,C) 

CALL SOLVE(MT.IPS,Z,BP.CP) 

DO 220 L=1,MT 

W'RITE(7,*) L.L*DPHI,CABS(C(L) + CP(L)),ATAN2(REAL(C(L) + CP(L)), 
* AIMAG(C(L) + CP(L)))/RAD 

ET= ET + RW(L)*C(L) + RW(L)*CP(L) 

220 EP = EP + RW(L + MT)*C(L) + RW(L + MT)*CP(L) 

ENDIF 

ET=UC*ET 

EP=UC*EP 

C 

C E-THETA IS CO-POL; E-PHI IS CROSS-POL 
C 

ANG(I)=THETA 
ECV(I) = CABS(ET) 

EXV(I) = CABS(EP) 

ECR=REAL(ET) 

ECI = AIM.^G(ET) 

EXR=REAL(EP) 

EXI = AIMAG(EP) 

PHC(I) = ATAN2(ECLECR+ l.E-20)/RAD 
PHX(I) = ATAN2(EXLEXR+ I.E-20)/RAD 
ECX = AMAXI(ECX,ECV(I).EXV(I)) 

500 CONTINUE 
DO 600 I=1,IT 

ECV(I) = AMAXI(ECV(I).1.E-10) 

EXV(I) = AMAX1(EXV(I),1.E-10) 

ECP(I) = (ECV(I)/ECX)**2 
EDP(I) = (EXV(I)/ECX)**2 
ECP(I) = AMAXI(ECP(I).. 00001) 

EDP(I) = AMAX1(EDP(I), .00001) 

ECP(I) = 1 0. *ALOG 1 0(ECP(I)) 

EDP(I)=10.*ALOG10(EDP(I)) 

600 CONTINUE 

SIGM A = (ECX**2)*CABS(UC)/(2. *BK) 

SIGDB = I0.*ALOGI0(SIGMA) 

WRITE(6,*) ’SIGMA, IN DB = ’,SIGMA,SIGDB 
DO 9000 L=I,IT 
WRITE(8,*) ANG(L),ECV(L) 

9000 CONTINUE 
900 CONTINUE 
STOP 
END 
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SUBROUTINE SOLVE(N,IPS,UL,B,X) 



SUBROUTINE TO SOLVE SYSTEM OF EQUATIONS WITH COMPLEX 
OEFFICIENTS. 

CALL ’DECOMP’ FIRST. (FROM MAUTZ AND HARRINGTON) 

COMPLEX UL(50000),B(500),X(500),SUM 
DIMENSION IPS(SOO) 

NP1=N + 1 
IP = IPS(1) 

X(1) = B(IP) 

DO 2 1 = 2, N 
IP = IPS(I) 

IPB = IP 
IM1=I-1 
SUM = (0.,0.) 

DO 1 J=1,IM1 

SUM = SUM + UL(IP)+X(J) 

1 IP=IP + N 

2 X(I) = B(IPB)-SUM 
K2 = N*(N-1) 

IP=IPS(N) + K2 
X(N) = X(N)/UL(IP) 

DO 4 IBACK = 2,N 
I = NP1-IBACK 

K2 = K2-N 
IPI=IPS(I) + K2 
IP1=I + 1 
SUM=(0.,0.) 

IP=IPI 

DO 3 J = IP1,N 
IP = IP + N 

3 SUM = SUM + UL(IP)*X(J) 

4 X(I) = (X(I)-SUM)/UL(IPI) 

RETURN 

END 

SUBROUTINE DECOMP(N,IPS,UL) 

SUBROUTINE TO DECOMPOSE SYSTEM OF EQUATIONS. 

FROM MAUTZ AND HARRINGTON. 



COMPLEX UL(50000), PIVOT, EM 
DIMENSION SCL(500),IPS(500) 

DO 5 I=1,N 
IPS(I) = I 
RN = 0. 

J1=I 

DO 2 J = 1,N 

ULM = ABS(REAL(UL(J 1))) + ABS(A1M AG(UL(J 1 ))) 
J1=J1 + N 
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IF(RN-ULM) 1,2,2 
1 RN=ULM 
2 CONTINUE 
SCL(I)=1./RN 
5 CONTINUE 
NM1 = N-1 
K2 = 0 

DO 17 K=1,NM1 
BIG=0. 

DO 11 I = K,N 
IP=IPS(I) 

IPK=IP+K2 

S1ZE=(ABS(REAL(UL(IPK))) + ABS(A1MAG(UL(IPK))))*SCL(IP) 

IF(SIZE-BIG) 11,11,10 

10 B1G = SIZE 
IPV = 1 

11 CONTINUE 
IF(IPV-K) 14,15,14 

14 J = IPS(K) 

1PS(K) = IPS(IPV) 

IPS(1PV) = J 

15 KPP = IPS(K) + K2 
PIVOT=UL(KPP) 

KP1=K+1 

DO 16 I = KP1,N 
KP = KPP 
IP = 1PS(1) + K2 
EM = -UL(1P)/P1V0T 
18 UL(IP) = -EM 
DO 16 J = KP1,N 
IP = 1P + N 
KP=KP+N 

UL(IP) = UL(IP) + EM *UL(KP) 

16 CONTINUE 
K2=K2 + N 

17 CONTINUE 
RETURN 
END 

SUBROUTINEZMATWW(NWIRES,NWl,NW2,RB.XH,YH,ZH,NT,XT.AT,AK,ZZ) 

IMPEDANCE ELEMENTS FOR LINEAR BASIS FUNCTIONS. 

COMPLEX CEXP,Z(200),ZZ(I5000),CON,CMPLX,EXP 

COMPLEX U0,SUM,SUM1,SUM2,SUM3,SUM4 

DIMENSION ZH(200),XT(128),AT(I28).XH(200),YH(200),UU(200) 

DIM ENSI ON D 1 (200) , S I (200) , C 1 (200) , ZS I (200) 

DIMENSION XS1(200),YS 1(200) 

DIMENSION CU(200),SU(200) 

INTEGER NT,NWIRES,NW1(4).NW2(4),NS(4) 

PI = 3. 1415926 
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PI2 = 2.*PI 
ETA = 377. 

BK = PI2 
U0 = (0.,0.) 

CON = (0. , 1 . )*BK*ETA/(4. *PI) 
A=AK 



DEFINE GEOMETRY TERMS FOR THE WIRE. XH,YH,ZH ARE ALL KNOWNS. 

DO 5 L=1,NWIRES 
NS(L) = NW2(L)-NW1(L)+ 1 
NSl =NW2(NW1RES)-NW1(1) 

NPS = NS1 + 1 
NTRlA = NPS-2 
DO 10 N = 2,NPS 
N0=N-1 
I = NW1(1) + N-1 
12 = 1-1 



AVERAGE VALUES 

ZS1(N0)= .5*(ZH(1) -f-ZH(I2)) 
XS1(N0) = .5*(XH(1) -H XH(I2}) 
YS1(N0) = .5*(YH(I) -h YH(I2)) 

DX = XH(I)-XH(I2) 

DY = YH(I)-YH(I2) 

Dl(N0) = SQRT(DX**2-i-DY**2) 
UU(N0) = ATAN2(DY,DX-(- l.E-5) 
CU(N0) = COS(UU(N0)) 

SU(N0) = SIN(UU(N0)) 

S1(N0) = DR/D1(N0) 

C1(N0) = DZ/D1(N0) 

10 CONTINUE 
IP=1 

WRITE(6,*) ’IP=MP 
DO 600 JQ=LNTRIA 

DOING II 



I = IP 
J=JQ 

CC = COS(UU(I)-UU(J)) 
TIP=1./D1(I) 
TJP=1./D1(J) 
Tl=Dl(I)/2. 

T2 = Dl(J)/2. 

SUM = U0 
DO 100 K=1,NT 
T=T1*XT(K) 

T1 = .5+T/D1(1) 
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X1 = XS1(I) + T*CU(I) 

Y1 = YS1(I) + T*SU(1) 

ZI = ZS1(I) 

DO 100 L=1,NT 
TP=T2*XT(L) 

TJ = .5 + TP/D1(J) 

XJ = XS1(J) + TP*CU(J) 

YJ = YS1(J) + TP*SU(J) 

ZJ = ZS1(J) 

RP = SQRT((XI-XJ)**2 + (YI-YJ)**2 + A**2) 

EXP=CEXP(CMPLX(0..-RP))/RP 

SUM = SUM + AT(L)*AT(K)*EXP*(TI‘n'J*CC-TIP*TJP) 

100 CONTINUE 
SUMl=SUM*Tl»CON*T2 

C 

C DOING 12 
C 

J=JQ+ 1 

CC = COS(UU(l)-UU(J)) 

T1P=1./D1(1) 

TJP=-1./D1(J) 

Tl=Dl(I)/2. 

T2 = Dl(J)/2. 

SUM=U0 
DO 101 K=1.NT 
T = T1*XT(K) 

TI = .5 + T/D1(I) 

XI = XS1(I) + T*CU(1) 

Y1 = YS1(I) + T*SU(1) 

ZI = ZS1(I) 

DO 101 L=1,NT 
TP = T2'^XT(L) 

TJ = .5-TP/D1(J) 

XJ = XS1(J)+TP*CU(J) 

YJ = YS1(J) + TP*SU(J) 

Z.1 = ZS1(J) 

RP = SQRT((X1-XJ)**2 + (Yl- YJ)**2 + A**2) 

EXP = CEXP(CMPLX(0. ,-RP))/RP 

SUM = SUM + AT(L)*AT(K)*EXP*(T1*TJ '^CC-T1P*TJP) 

101 CONTINUE 
SUM2=SUM-Tl*CON-n'2 

C 

C DOING 13 
C 

I = IP+1 
J=JQ 

CC = COS(UU(I)-UU(J)) 

TIP = -1./D1(1) 

TJP=1./D1(J) 

Tl=Dl(I)/2. 
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T2 = Dl(J)/2. 

SUM = U0 

DO 102 K=1,NT 
T=T1*XT(K) 

TI = .5-T/D1(I) 

XI = XS1(I) + T*CU(I) 

YI = YSl(I)-fT*SU(I) 

ZI = ZS1(I) 

DO 102 L=1,NT 
TP=T2*XT(L) 

TJ = .5+TP/D1(J) 

XJ = XSl(J)-fTP‘"CU(J) 

YJ = YS1(J) + TP*SU(J) 

ZJ = ZS1(J) 

RP = SQRT((XI-XJ)**2 + ( YI- YJ)**2 + A**^2) 

EXP = CEXP(CM PLX(0. ,-RP))/RP 

SUM = SUM + AT(L)*AT(K)*EXP*(TI*TJ*CC-TIP^JP) 

102 CONTINUE 

SUM3 = SUM *T1 ’"CON*T2 
C 

C DOING 14 
C 

J=JQ-hl 

CC = COS(UU(I)-UU(J)) 

TIP = -1./D1(I) 

TJP = -1./D1(J) 

Tl=Dl(I)/2. 

T2 = Dl(J)/2. 

SUM = U0 

DO 103 K=1,NT 
T = T1*XT(K) 

TI = .5-T/D1(1) 

XT = XS1(I) + T*CU(I) 

YI = YS1(I) + T*SU(I) 

ZI = ZS1(I) 

DO 103 L=1,NT 
TP=T2*XT(L) 

TJ = .5-TP/D1(J) 

XJ = XS1(J) + TP*CU(J) 

YJ = YS1(J) + TP*SU(J) 

ZJ = ZS1(J) 

RP = SQRT((XI-XJ)*’^2 + (YI-YJ)**2 + A**2) 

EXP = CEXP(CMPLX(0.,-RP))/RP 

SUM = SUM + AT(L)*AT(K)*EXP*(TI *TJ*CC-TIP*TJP) 

103 CONTINUE 

SUM4 = SUM*T1*C0N^2 
C 

C IMPEDANCE ELEMENT FOR IP,JQ 
C 

KK=(JQ-1)*NTRIA + IP 
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Z(JQ) = (SUM 1 + SUM2 + SUM 3 + SUM4) 

600 CONTINUE 
Z(NTRIA)=Z(2) 

C 

C FILL THE ENTIRE Z MATRIX USING SYMMETRY PROPERTIES 
C ROW INDEX. I; COL INDEX, J 
C 

DO 12 I=1,NTRIA 
DO 12 J = 1,NTR1A 
K = (M)*NTRIA+J 
ZZ(K) = U0 
U = IABS(I-J) 

IF(U.GT.NTRIA) GO TO 12 

U1=U+1 

ZZ(K) = Z(U1) 

12 CONTINUE 
CLOSE(2) 

RETURN 

END 

SUBROUTINE PLANEW(NP.XH.YH,ZH.THR,PHR,R) 

C 

C PLANE WAVE EXCITATION VECTOR ELEMENTS FOR WIRE AND 
C INCIDENCE DIRECTION IS (THR.PHR). 

C WIRE LIES IN THE X-Y PLANE (Z=0) 

C 

COMPLEX U0,C,R(2000),CEXP,EXP,FI1,F12,S1.DLCMPLX 
DIMENSION ZH(500),XH(500),YH(500) 

MP2 = NP-1 
MT2 = NP-2 
U0=(0.,0.) 

CC = COS(THR) 

SS = S1N(THR) 

CP = COS(PHR) 

SP=SIN(PHR) 

UP = SS*CP 
VP=SS*SP 
DO 12 IP=I,MP2 
II = IP 
1 = 11 + 1 

ZS = .5*(ZH(I) + ZH(II)) 

XS = .5*(XH(I) + XH(I1)) 

YS = .5*(YH(I) + YH(1I)) 

DX = XH(I)-XH(II) 

DY = YH(1)-YH(II) 

D 1 = SQRT(DX**2 + D Y**2) 

SU = DY/D1 
CU = DX7D1 

C FOR WIRES IN THE XY PLANE SIN(V)= 1 AND COS(V) = 0 
SV=1.0 
CV = 0.0 
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WIRE SEGMENT CALCULATIONS 

A = UP*CU + VP*SU 
B = UP*XS + VP*YS 
C = CMPLX(0.,A) 
EXP=CEXP(CMPLX(0.,B)) 

AA = CC*(CU*CP + SU*SP) 

BB = SU*CP-SP*CU 
PSI = Dl*A/2. 

1F(PSI.NE.0.) GO TO 60 
SINC=1. 

GO TO 61 

60 S1NC = S1N(PSI)/PSI 

61 COSP=COS(PSI) 
FIl=SINC*Dl*EXP/2. 

F12 = (0.,0.) 

lF(ABS(A).LT.l.E-4) GO TO 62 
CSP = COSP-SlNC 
IF(ABS(CSP).LT.l.E-4) GO TO 62 
F12 = EXP/C*CSP 

62 CONTINUE 
SI = F11+F12 
D1 = FI1-F12 

R-WIRE-THETA 

1F(IP.EQ.MP2) GO TO 10 
R(1P) = AA*S1 
R(1P + MT2) = BB*S1 
10 CONTINUE 

R-WIRE-PHI 

14 lF(lP.EQ.l)GO TO 12 
R(IP-1) = R(1P-1) + AA*D1 
R(IP-1+MT2) = R(1P-1+MT2) + BB*D1 
12 CONTINUE 
210 RETURN 
END 
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C MAIN PROGRAM *HARLOOP.F* 

C PLANE WAVE SCATTERING FROM A CIRCULAR LOOP IN THE Z PLANE. 

C RADAR CROSS SECTION CALCULATION ***** 

C USING HARRINGTON’S FORMULATION FROM THE BOOK ’FIELD COMP. BY 
C MM’ (P.83 TO 95) 

C 

COMPLEX Z(5000),E(250),C(250),EX,EC,ET,EP,RW(1000),UC,EPHI(500) 
COMPLEX CPHI(500) 

DIMENSION ECP(500),IPS(250),ANG(500),EDP(500),XT{300),AT(300) 
DIMENSION ECV(500),EXV(500),PHC(500),PHX(500) 

DATA PI/3.14159265/ 

DATA IPRINT/1/,ITEST/1/ 

RAD = PI/1 80. 

ECX = 0. 

BK = 2.*PI 
ETA = 377. 

UC = (0,-1 )*ETA*BK/(4. *PI) 

PHIRD = 0. 

OPEN( 1 , FILE = ’ PAR AM LST. D AT’ ) 

READd,*) ANGLE, DT, START, STOP, A,B,SEG,AHARR,GCONST,XLOC,XlPOL 
IPOL=INT(XIPOL) 

LOC = INT(XLOC) 

CLOSE(l) 

NM = INT(AHARR) 

CON = (377.*BK)**2/2./BK 
OPEN(2,FILE = ’OUTGLEG’) 

READ(2,*) NT 
DO 1 I=1,NT 
READ(2,*) XT(I),AT(1) 

1 CONTINUE 

OPEN(8,FlLE = ’OUTHARR.DAT’) 

OPEN(7,FILE = ’IHARROUT.DAT') 

NROW = 2*NM-i-l 
WRITE(6,1300) B,A.NROW,NT 

1300 FORMAT(//,5X, ’PLANE WAVE SCATTERING BY A CIRCULAR LOOP’,/ 

* ,5X, ’USING METHOD IN HARRINGTONS MM TEXT BOOK’,/ 

* ,5X,’LOOP RADIUS (WAVL) = ’,F8.4,/,5X,’WIRE RADIUS (U^AVL) = ’, 

* F8.4,/,5X,’NUMBER OF AZIMUTHAL MODES (INCLUDING ZERO) = ’,I3, 

* /,5X, ’NUMBER OF INTEGRATION POINTS IN PHI=’,I4) 

COMPUTE IMPEDANCE MATRIX ELEMENTS 

WRITE(6,*) ’CALLING ZMAT’ 

CALL ZMATWW'(NM,A,B,NT,XT,AT,Z) 

WRITE(6,*) ’WIRE IMPEDANCE COMPUTED’ 

CALL DECOMP(NROW,IPS,Z) 

WRITE(6,*) ’Z DECOMPOSED’ 

BEGIN FIELD CALCULATIONS 
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PHR0 = PHIRD*RAD 
IT=INT((STOP>START)/DT) -f 1 
WRITE(7,‘»‘) IT,NROW 
DO 500 1=1, IT 

THETA = FLOAT(l- 1 )*DT -f START 
WRITE{6,*) THETA = \THETA 
THR=THETA*RAD 
PHR=PHR0 

IF(THETA.LT.180.) GO TO 99 
THR = (360.-THETA)*RAD 
PHR = PHR0-fPI 
99 CONTINUE 
ET = (0.,0.) 

EP = (0.,0.) 

CALL PLANEW(NM,B,THR,PHR,RW) 

C TRANSMIT VECTOR ELEMENTS ARE TRANSPOSED FORMS OF RECEIVE 
VECTOR 

C ALSO THE THETA COMPONENT GETS A NEGATIVE SIGN 
IF(LOC .EQ. 0) THEN 
IF(IPOL.E0-l)THEN 
C 

C THETA POLARIZED INCIDENT WAVE (IPOL=l) 

C 

DO 101 L=LNROW 

101 E(NROW-L-h 1) = RW(L) 

ELSE 

C 

C PHI POLARIZED INCIDENT WAVE (IPOL=2) 

C 

DO 102 L=LNROW 

102 E(NROW-L+ l) = RW(L + NROW) 

ENDIF 

WRITE(6,*) 'CALLING SOLVE' 

CALL SOLVE(NROW,IPS,Z,E,C) 

WRITE(6,*) 'RETURNED FROM SOLVE' 

DO 210 L=LNROW 
WTUTE(7,*) C(L) 

ET=ET + RW(L)*C(L) 

210 EP = EP-f R W(L -f NROW) *C(L) 

ELSE 

C 

C E-THETA IS CO-POL; E-PHI IS CROSS-POL 
C 

DO 221 L=LNROW 

22 1 E(NROW-L + 1 ) = RW(L) 

DO 222 L=l,NROW 

222 EPHI(NROW-L-f l) = RW(L + NROW)*CEXP(CMPLX(0.,PI/2.)) 

CALL SOLVE(NROW,IPS,Z,E,C) 

CALL SOLVE(NROW,IPS,Z,EPHI,CPHI) 

DO 220 L=l,NROW 
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WR1TE(7,*) C(L) + CPHI(L) 

ET= ET + R W(L)*C(L) + R W(L)*CPHI(L) 

220 EP = EP + RW(L + NROW)*C(L) + RW(L + NROW)*CPHl(L) 

ENDIF 
EC = UC*ET 
EX = UC*EP 
ANG(I)=THETA 
ECV(I) = CABS(EC) 

EXV(1) = CABS(EX) 

ECR = REAL(EC) 

ECI = A1MAG(EC) 

EXR = REAL(EX) 

EX1 = AIMAG(EX) 

PHC(1 ) = AT AN2(ECI , ECR + 1 . E-20)/R AD 
PHX(1) = ATAN2(EXI,EXR + 1 ,E-20)/RAD 
ECX = AMAX1(ECX,ECV(1),EXV(I)) 

500 CONTINUE 
DO 600 1=1, IT 

ECV(l) = AM AX 1 (EC V(l), 1 . E- 1 0) 

EXV(I) = AMAX1(EXV(I),1.E-10j 
ECP( l) = (ECV(])/ECX)**2 
EDP(I) = (EXV(I)/ECX)**2 
ECP(I) = AMAX1(ECP(1),. 00001) 

EDP(I) = AM AX 1(EDP(1), .00001 ) 

ECP(I)= 10.*ALOG10(ECP(1)) 

EDP(I)= 10.*ALOG10(EDP(I)) 

600 CONTINUE 

SIGMA = (ECX**2)*4.*PI 
SIGDB= 10.*ALOG10(S1GMA) 

WRITE(6,*) ’BACKSC.MTER, IN DB = ’, SIGMA, SIGDB 
OPEN(2,FILE= ’RCS.DAT’) 

WRITE(2,*) SIGMA, SIGDB 
CLOSE(2) 

C 

C PRINT FIELD POINTS 
C 

DO 9000 L=1,IT 
WRITE(8,*) ANG(L),ECV(L) 

5016 FORM AT(5X,F6. 1 ,3X,2(F8.4,3X,F7. 1 ,3X,F7.2,3X)) 

9000 CONTINUE 
900 CONTINUE 
9998 STOP 
END 

SUBROUTINE SOLVE(N,IPS,UL,B,X) 

COMPLEX UL(5000),B(250),X(250),SUM 
DIMENSION IPS(250) 

NP1=N + 1 
IP = IPS(1) 

X(1)=B(IP) 

DO 2 I=2,N 
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IP = IPS(I) 

IPB = IP 
IM1=I-1 
SUM = (0.,0.) 

DO 1 J = 1,IM1 

SUM = SUM + UL(IP)*X(J) 

1 IP=IP+N 

2 X(1) = B(IPB)-SUM 
K2=N*(N-1) 

IP=IPS(N) + K2 
X(N) = X(N)/UL(IP) 

DO 4 IBACK = 2,N 
I = NP1-IBACK 
K2=K2-N 
IPI=IPS(I) + K2 
IP1=I + 1 

SUM = (0.,0.) 

1P=IPI 

DO 3 J = IP1,N 
IP = IP + N 

3 SUM = SUM + UL(IP)*X(J) 

4 X(1) = (X(I)-SUM), UL(IP1) 

RETURN 

END 

SUBROUTINE DECOMP(N,IPS,UL) 

COMPLEX UL(5000), PIVOT, EM 
DIMENSION SCL(250),IPS(250) 

DO 5 1=1, N 
IFS(I) = I 
RN=0. 

J1=I 

DO 2 J=1,N 

ULM = ABS(REAL(UL(J 1))) + ABS(AIM AG(UL(J 1 ))) 

J1=J1+N 
IF(RN-ULM) 1,2,2 
1 RN=ULM 
2 CONTINUE 
SCL(I)=1./RN 

5 CONTINUE 
NM1 = N-1 
K2 = 0 

DO 17 K=1,NM1 
BIG = 0. 

DO 11 I = K,N 
IP = IPS(1) 

IPK = IP + K2 

SIZE = (ABS(REAL(UL(IPK))) + ABS(AIM AG(UL(IPK))))*SCL(IP) 
IF(SIZE-BIG) 11,11,10 
10 BIG = SIZE 
IPV = I 
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11 CONTINUE 
IF(IPV-K) 14,15,14 

14 J = IPS(K) 

IPS(K) = IPS(IPV) 

1PS(IPV)=J 

15 KPP = IPS(K) + K2 
PIVOT = UL(KPP) 

KP1=K+1 

DO 16 I=KP1,N 
KP=KPP 
IP=IPS(I)+K2 
EM=-UL(IP)/P1V0T 
18 UL(IP) = -EM 
DO 16 J=KP1,N 
IP=IP + N 
KP=KP+N 

UL(IP) = UL(IP) + EM *UL(KP) 

16 CONTINUE 
K2=K2 + N 

17 CONTINUE 
RETURN 
END 

SUBROUTINE ZM ATVAV'(NM,A,B,NT,XT,AT,Z) 

C 

C *** MODS FOR LOOP - USING HARRINGTON’S TEXT BOOK EQUATIONS AS A 
C CHECK FOR OTHER METHODS. NM IS THE NUMBER OF AZIMUTHAL MODES 
USED 
C 

COMPLEX Z(5000).CON.CMPLX,FK 
DIMENSION XT(300),AT(300) 

PI = 3. 1415926 
BK = 2.*PI 

CON = CMPLX(0. ,PI*377. *BK*B) 

NR0W = 2.*NM + 1 
DO 10 I = -NM,NM 
DO 10 J=-NM,NM 
U=J + NM + 1 + (KNM)*NRO\V 
Z(IJ)=(0.,0.) 

10 CONTINUE 
C 

C ONLY DIAGONAL ELEMENTS ARE NONZERO. ALTHOUGH SYMMETRY EXISTS 
BETWEEN 

C Z(-N,-N) AND Z(N,N) IT IS NOT BEING USED. 

C 

DO 20 I = -NM,NM 
J=I + NM + l+(I + NM)*NROW 
IP=I + 1 
IM = I-1 

Z(J) = (.5*FK(IM.B,A,NT,XT,AT) + .5*FK(IP,B,A,NT,XT,AT)- 
* (I/BK/B)**2*FK(I,B,A,NT,XT,AT))*C0N 
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20 CONTINUE 
RETURN 
END 

SUBROUTINE PLANEW(N,B,THR,PHR,R) 

C 

C PLANE WAVE RECEIVE VECTOR ELEMENTS FOR WIRE USING THE 
C FORMULATION FROM HARRINGTON’S BOOK. N IS THE NUMBER OF 
C AZIMUTHAL MODES. NOTE THAT B(N) = R(-N) (B IS EXCITATION AND 
C R IS RECEIVE). 

C 

COMPLEX R(I000),CEXP,EXP,CMPLX 
PI = 3. 14159 
BK = 2.*PI 
CT=COS(THR) 

ST=SIN(THR) 

RR = BK*B*ST 
NROW = 2*N+l 

C DO THETA RECEIVE COMPONENTS FIRST 
DO 10 I = -N,N 
IP = I + 1 
IM = M 

EXP = CEXP(CMPLX(0.,I*PHR)) 

Ra + N+1) = -PPB*(0.,1.)**I*EXP*(BESSJ(IP,RR) + BESSJ(IM,RR))*CT 
10 CONTINUE 

C NOW DO PHI RECEIVE COMPONENTS 
DO 20 I = -N,N 
IP = I+1 
IM = M 

EXP = CEXP(CMPLX(0.,I*PHR)) 

R(I N -4- 1 + NROW) = PI*B*(0. , 1 . )**(! + 1 )*EXP*(BESSJ(IP,RR) 
-BESSJ(IM,RR)) 

20 CONTINUE 
RETURN 
END 

COMPLEX FUNCTION FK(N,B,A,NT,XT,AT) 

COMPLEX SUM,EXP1,EXP2,CEXP,CMPLX 
DIMENSION XT(300),AT(300) 

PI = 3. 14159 
BK = 2.’*'PI 
CC=1,/BK 
Pl = (2.*^PL0.)/2. 

P2 = (2.^PI*fO.)/2. 

SUM = (0.,0.) 

DO 10 1=1, NT 
PHI = P1*XT(I) + P2 

RR = SQRT(4. *SIN(PHI/2. )**2 + (A/B)*^^2) 

EXPl =CEXP(CMPLX(0.,-BK*B^RR)) 

EXP2 = CEXP(CMPLX(0. ,-N*PHI)) 

SUM = SUM + AT(I)*EXP 1 *EXP2/RR 
10 CONTINUE 
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FK = SUM*P1*CC 

RETURN 

END 

FUNCTION BESSJ(NN,X) 

C RETURNS THE BESSEL FUNCTION B OF ORDER N (> 1) AND REAL 
C ARGUMENT X. 

PARAMETER (lACC = 40.B1GNO = 1 .E lO.BIGNI = 1 . E- 10) 

IF(NN.LT.O) N = -NN 
IF(NN.GE.O) N=NN 
KC = 3 

IF(N.EQ.O) KC = 1 
IF(N.EQ.l) KC=2 
GO TO (1,2,3), KC 

1 BESSJ=BESSJO(X) 

GO TO 4 

2 BESSJ = BESSJ1(X) 

GO TO 4 

3 BESSJ = 0. 

IF(ABS(X).LT.l.E-5) GO TO 4 
TOX = 2./X 

IF(X.GT.FLOAT(N)) THEN 
BJM = BESSJ0(X) 

BJ=BESSJ1(X) 

DO 11 J = 1,N-1 
BJP = J*TOX*BJ-BJM 
BJM=BJ 
BJ = BJP 

1 1 CONTINUE 
BESSJ = BJ 

ELSE 

M =2*((N + INT(SQRT(FLOAT(lACC-"N))))/2) 

BESSJ = 0. 

JSUM = 0. 

SUM=0. 

BJP = 0. 

BJ = 1. 

DO 12 J = M,1,-1 
BJM=J*TOX*BJ-BJP 
BJP = BJ 
BJ = BJM 

IF(ABS(BJ).GT.BlGNO) THEN 
BJ = BJ*BIGNI 
BJP = BJP*BIGNI 
BESSJ = BESSJ*B1GNI 
SUM = SUM*BIGNI 
ENDIF 

IF(JSUM.NE.O) SUM = SUM + BJ 
JSUM = 1-JSUM 
IF(J.EQ.N) BESSJ = BJP 

12 CONTINUE 
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SUM = 2.*SUM-BJ 
BESSJ = BESSJ/SUM 
ENDIF 

4 CONTINUE 

IF(NN.LT.O) BESSJ = (-1.)**N*BESSJ 
RETURN 
END 

FUNCTION BESSJO(X) 

BESSEL FUNCTION OF 0 ORDER, REAL ARGUMENT X 
(SEE ’NUMERICAL RECIPES’. P.172) 

REALMS Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 

* S1,S2,S3,S4,S5,S6 

DATA Pl,P2,P3,P4,P5/l.D0,-.109862827D-2,.2734510407D-4, 

* -.2073370639D-5, . 209388721 lD-6/ 

DATA Ql,Q2,Q3,Q4,Q5/-.1562499995D-l,.1430488765D-3, 

* -.6911147651D-5,.7621095161D-6,-.934945152D-7/ 

DATA R1,R2,R3,R4,R5,R6/57568490574.D0, -13362590354. DO, 

* 651619640.7D0.-1 1214424. 18D0,77392.33017D0,-184.9052456D0/ 
DATA S1,S2,S3,S4,S5,S6/57568490411.D0,1029532985.D0, 

* 9494680.7 1 8D0,59272.64853D0,267 .85327 12D0, 1 .DO/ 
IF(ABS(X).LT.8.)THEN 

Y = X**2 

BESSJ0 = (R 1 -t- Y*(R2 -H Y»(R3 -i- Y*(R4 -I- Y*(R5 -i- Y*R6)))))/ 

* (SI Y*(S2 + Y*{S3 -I- Y*(S4 -I- Y*(S5 -i- Y’^S6))))) 

ELSE 

AX = ABS(X) 

Z = 8./AX 

Y = Z**2 

XX = AX-.785398164 

BESSJO = SQRT(.636619772/AX)*(COS(XX)*(Pl+Y*(P2-i-Y*(P3-i- 

* Y*(P4 -I- Y*P5))))-Z*SIN(XX)*(Q 1 -I- Y*(Q2 -i- Y*(Q3 -l- 

* Y*(Q4-l-Y*Q5))))) 

ENDIF 
RETURN 
END 

FUNCTION BESSJl(X) 

C 

C BESSEL FUNCTION B OF ORDER 1, REAL ARGUMENT X 
C (SEE 'NUMERICAL RECIPES’, P.173) 

C 

REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 

* S1,S2,S3,S4,S5,S6 

DATA Pl,P2,P3,P4,P5/l.D0,.183105D-2,-.3516396496D-4, 

* .2457520174D-5.-.20337019D-6/ 

DATA Ql,Q2,Q3,Q4,Q5/.04687499995D0,-.2002690873D-3, 

* .8449199096D-5.-.99228987D-6,. 105787412D-6/ 

DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0, 

* 242396853. 1 D0.-297261 1 .439D0, 15704.4826D0.-30. 16036606D0/ 
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DATA SI, S2,S3.S4,S5.S6/144725228442.D0, 2300535 178. DO, 

* 1 8583304. 74D0,99447.43394D0,376.999 1397D0, 1 .DO/ 
IF(ABS(X).LT.8.) THEN 

Y = X**2 

BESSJ 1 = X*(R 1 + Y*(R2 + Y*(R3 + Y*(R4 + Y*(R5 + Y*R6)))))/ 

* (S 1 + Y*(S2 + Y*(S3 + Y*(S4 + Y*(S5 + Y*S6))))) 

ELSE 

AX = ABS(X) 

Z = 8./AX 

Y = Z**2 

XX = AX-2.356194491 

BESSJ 1 = SQRT(.6366 19772/AX)*(COS(XX)*(P 1 + Y*(P2 + Y*(P3 + 

* Y*(P4 + Y*P5))))-Z*SIN(XX)*(Q1 + Y*(Q2 + Y*(Q3 + 

* Y*(Q4 + Y*Q5)))))*SIGN(1.,X) 

ENDIF 

RETURN 

END 
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